﻿#include "Algorithms.hpp"
#include "Proc/OiFitLine.hpp"

#define CGAL_NDEBUG
#define CGAL_NO_AUTOLINK
#define CGAL_CFG_NO_NEXTAFTER

#define BOOST_ALL_NO_LIB
#define BOOST_ERROR_CODE_HEADER_ONLY

#undef CGAL_error
#define CGAL_error()

#include <CGAL/assertions.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/intersections.h>
#include <CGAL/Exact_circular_kernel_2.h>
#include <CGAL/Circular_kernel_intersections.h>
#include <CGAL/iterator.h>
#include <CGAL/Polygon_2.h>


typedef CGAL::Simple_cartesian<double> K;
typedef CGAL::Exact_circular_kernel_2 C;
typedef std::pair<CGAL::Circular_arc_point_2<C>, unsigned> AP;

namespace Oi
{
    namespace PLScannerUtil
    {
        namespace
        {
            inline K::Point_2& as(QPointF& point)
            {
                return *reinterpret_cast<K::Point_2*>( &point );
            }
            inline C::Point_2& asC(QPointF& point)
            {
                return *reinterpret_cast<C::Point_2*>( &point );
            }
            inline QPointF& as(K::Point_2& point)
            {
                return *reinterpret_cast<QPointF*>( &point );
            }
            inline QPointF to(const K::Point_2& point)
            {
                return QPointF(point.x(), point.y());
            }
            inline double euler(const QPointF& dir)
            {
                return sqrt(dir.x() * dir.x() + dir.y() * dir.y());
            }
        }

        inline double pointInLine(const QPointF& p, const QLineF& l)
        {
            auto p1 = l.p1();
            auto p2 = l.p2();
            return ( p2.y() - p1.y() ) * p.x() + ( p1.x() - p2.x() ) * p.y() + p2.x() * p1.y() - p1.x() * p2.y();
        }

        double CalcDistance(const QPointF& p1, const QPointF& p2)
        {
            return sqrt(pow(p1.y() - p2.y(), 2) + pow(p1.x() - p2.x(), 2));
        }

        double CalcDistance(const QPointF& p, const QLineF& l)
        {
            auto foot = perpendicularFoot(p, l);
            if (foot.isNull())
            {
                return 0.0;
            }
            return CalcDistance(p, foot);
        }

        QPointF perpendicularFoot(const QPointF& p, const QLineF & l)
        {
            if(!qFuzzyIsNull(l.length()))
            {
                auto dotmultiply = ( p.x() - l.p1().x() ) * ( l.p2().x() - l.p1().x() ) + ( p.y() - l.p1().y() ) * ( l.p2().y() - l.p1().y() );
                auto relation = dotmultiply / ( ( l.p2().x() - l.p1().x() ) * ( l.p2().x() - l.p1().x() ) + ( l.p2().y() - l.p1().y() ) * ( l.p2().y() - l.p1().y() ) );

                return QPointF(l.p1().x() + relation * ( l.p2().x() - l.p1().x() ), l.p1().y() + relation * ( l.p2().y() - l.p1().y() ));
            }
            else
            {
                return QPointF();
            }
        }

        double judgeSide(const QPointF& p1, const QPointF& p2, const QLineF & l)
        {
            return pointInLine(p1, l) * pointInLine(p2, l);
        }

        bool DFT_1D(const cv::Mat line, cv::Mat_<float>& amp, cv::Mat_<float>& phase)
        {
            if(line.channels() != 1 || (line.cols != 1 && line.rows != 1))
            {
                return false;
            }

            cv::Mat dftMats[] = { line, cv::Mat::zeros(line.size(), CV_32F) };
            cv::Mat dftMat;
            cv::merge(dftMats, 2, dftMat);
            cv::dft(dftMat, dftMat);

            std::vector<cv::Mat> channels;
            cv::split(dftMat, channels);

            cv::Mat_<float> mRe = channels[0];
            cv::Mat_<float> mIm = channels[1];

            cv::magnitude(mRe, mIm, amp);

            amp /= amp.cols / 2;

            phase = cv::Mat_<float>(amp.size());

            float* mReD = (float*)mRe.data;
            float* mImD = (float*)mIm.data;
            float* phaseD = (float*)phase.data;

            for(int i = 0; i < phase.cols; i++, mReD++, mImD++, phaseD++)
            {
                if(*mReD == 0)
                {
                    if(*mImD == 0)
                    {
                        *phaseD = 0;
                    }
                    else
                    {
                        if(*mImD > 0)
                        {
                            *phaseD = M_PI / 2.0;
                        }
                        else
                        {
                            *phaseD = -M_PI / 2.0;
                        }
                    }
                }
                else
                {
                    *phaseD = atan2f(*mImD, *mReD);
               
                }
            }

            return true;
            
        }
        template<typename T>
        void  _interpolation( cv::Mat src1, cv::Mat index, cv::Mat& dst, int direct)
        {
            //index恒为竖向的
            cv::Mat src;
            src1.copyTo(src);
            switch(direct)
            {
            case 0://src和dst竖向
                if (src.cols != 1||index.cols!=1)
                {
                    return;
                }
                if (src.rows == index.rows)
                {
                    std::list<int> interps;
                    std::list<int> poss;
                    std::list<float> values;

                    for (int i = 0; i < src.rows; i++)
                    {
                        if (index.ptr<T>(i)[0] == 0)
                        {
                        }
                        else
                        {
                            for (int j = 0; j < i; j++)
                            {
                                src.ptr<T>(j)[0] = src.ptr<T>(i)[0];
                                index.ptr<T>(j)[0] = 1;
                            }
                            break;
                        }
                    }

                    for (int i = src.rows - 1; i >= 0; i--)
                    {
                        if (index.ptr<T>(i)[0] == 0)
                        {
                        }
                        else
                        {
                            for (int j = src.rows - 1; j > i; j--)
                            {
                                src.ptr<T>(j)[0] = src.ptr<T>(i)[0];
                                index.ptr<T>(j)[0] = 1;
                            }
                            break;
                        }
                    }

                    for (int i = 0; i < src.rows; i++)
                    {
                        if (index.ptr<T>(i)[0] == 0)
                        {
                            interps.push_back(i);
                        }
                        else
                        {
                            if (poss.empty() || poss.back() + 1 < i)
                            {
                                poss.push_back(i);
                                values.push_back(src.ptr<T>(i)[0]);
                            }
                            else
                            {
                                poss.back() = i;
                                values.back() = src.ptr<T>(i)[0];
                            }

                            if (!interps.empty() && poss.size() == 2)
                            {
                                auto len = poss.back() - poss.front();
                                for (auto item : interps)
                                {
                                    auto len0 = item - poss.front();
                                    auto len1 = item - poss.back();
                                    src.ptr<T>(item)[0] = (values.back() * len0 - values.front() * len1) / len;
                                }
                                poss.pop_front();
                                values.pop_front();
                                interps.clear();
                            }
                        }
                    }
                    dst = src.clone();
                }
                else {
                    return;
                }
                
                break;
            case 1://src和dst横向
                if (src.rows != 1 || index.cols != 1)
                {
                    return;
                }
                if (src.cols == index.cols)
                {
                    std::list<int> interps;
                    std::list<int> poss;
                    std::list<float> values;
                    T* sp = (T*)src.ptr(0);
                    for (int i = 0; i < src.cols; i++)
                    {
                        if (index.ptr<T>(i)[0] == 0)
                        {
                        }
                        else
                        {
                            for (int j = 0; j < i; j++)
                            {
                                sp[j] = sp[i];
                                index.ptr<T>(j)[0] = 1;
                            }
                            break;
                        }
                    }

                    for (int i = src.cols - 1; i >= 0; i--)
                    {
                        if (index.ptr<T>(i)[0] == 0)
                        {
                        }
                        else
                        {
                            for (int j = src.cols - 1; j > i; j--)
                            {
                                sp[j] = sp[i];
                                index.ptr<T>(j)[0] = 1;
                            }
                            break;
                        }
                    }

                    for (int i = 0; i < src.cols; i++)
                    {
                        if (index.ptr<T>(i)[0] == 0)
                        {
                            interps.push_back(i);
                        }
                        else
                        {
                            if (poss.empty() || poss.back() + 1 < i)
                            {
                                poss.push_back(i);
                                values.push_back(sp[i]);
                            }
                            else
                            {
                                poss.back() = i;
                                values.back() = sp[i];
                            }

                            if (!interps.empty() && poss.size() == 2)
                            {
                                auto len = poss.back() - poss.front();
                                for (auto item : interps)
                                {
                                    auto len0 = item - poss.front();
                                    auto len1 = item - poss.back();
                                    sp[item] = (values.back() * len0 - values.front() * len1) / len;
                                }
                                poss.pop_front();
                                values.pop_front();
                                interps.clear();
                            }
                        }
                    }
                    dst = src.clone();
                }
                break;
            case 2://src和dst横竖
                break;
            case 3://src和dst竖横
                break;
            default:
                qDebug() << "所选方式超出范围！";
                break;
            }
        }
        void OI_PLSCANNER_UTIL_API interpolation(const cv::Mat src, cv::Mat index, cv::Mat& dst, int direct)
        {
            if (src.channels() != 1)
            {
                return;
            }
            if (src.depth() == CV_8U)
            {
                _interpolation<uchar>(src,index, dst, direct);
            }
            else if (src.depth() == CV_32F)
            {
                _interpolation<float>(src,index, dst, direct);
            }
        }
        template<typename T>
        void _reduceAvg(const cv::Mat src, cv::Mat_<float>& dst, int direct, bool fillBlock)
        {
            auto dd = (float*)dst.data;
            if(direct == 0)
            {
                std::vector<float> rSum;
                std::vector<int> rCount;
                rSum.resize(src.cols);
                rCount.resize(src.cols);

                for(int j = 0; j < src.cols; j++)
                {
                    rSum[j] = 0.0;
                    rCount[j] = 0;
                }

                for(int i = 0; i < src.rows; i++)
                {
                    T* sp = (T*)src.ptr(i);
                    
                    for(int j = 0; j < src.cols; j++)
                    {
                        if(sp[j] > -1000.0 && sp[j] < +1000.0)
                        {
                            rSum[j] += sp[j];
                            rCount[j]++;
                        }
                    }
                }

                for(int j = 0; j < src.cols; j++)
                {
                    dd[j] = rSum[j] / rCount[j];
                }
            }
            else 
            {
                for(int i = 0; i < src.rows; i++)
                {
                    T* sp = (T*)src.ptr(i);

                    double sum = 0.0;
                    int count = 0;

                    for(int j = 0; j < src.cols; j++)
                    {
                        if(sp[j] > -1000.0 && sp[j] < +1000.0)
                        {
                            sum += sp[j];
                            count++;
                        }
                    }

                    dd[i] = sum / count;
                }
            }

            if (fillBlock)
            {
                std::list<int> interps;
                std::list<int> poss;
                std::list<float> values;

                for(int i = 0; i < dst.size[0]; i++)
                {
                    if(isnan(dd[i]) || dd[i] < -1000.0 || dd[i]  > 1000.0)
                    {
                    }
                    else
                    {
                        for(int j = 0; j < i; j++)
                        {
                            dd[j] = dd[i];
                        }
                        break;
                    }
                }

                for(int i = dst.size[0] - 1; i >= 0; i--)
                {
                    if(isnan(dd[i]) || dd[i] < -1000.0 || dd[i]  > 1000.0)
                    {
                    }
                    else
                    {
                        for(int j = dst.size[0] - 1; j > i; j--)
                        {
                            dd[j] = dd[i];
                        }
                        break;
                    }
                }

                for(int i = 0; i < dst.size[0]; i++)
                {
                    if(isnan(dd[i]) || dd[i] < -1000.0 || dd[i]  > 1000.0)
                    {
                        interps.push_back(i);
                    }
                    else
                    {
                        if(poss.empty() || poss.back() + 1 < i)
                        {
                            poss.push_back(i);
                            values.push_back(dd[i]);
                        }
                        else
                        {
                            poss.back() = i;
                            values.back() = dd[i];
                        }

                        if(!interps.empty() && poss.size() == 2)
                        {
                            auto len = poss.back() - poss.front();
                            for(auto item : interps)
                            {
                                auto len0 = item - poss.front();
                                auto len1 = item - poss.back();
                                dd[item] = ( values.back() * len0 - values.front() * len1 ) / len;
                            }
                            poss.pop_front();
                            values.pop_front();
                            interps.clear();
                        }
                    }
                }
            }
            else
            {
                for(int i = 0; i < dst.size[0]; i++)
                {
                    if(isnan(dd[i]) || dd[i] < -1000.0 || dd[i]  > 1000.0)
                    {
                        dd[i] = -1000.0;
                    }
                }
            }
        }

        void reduceAvg(const cv::Mat src, cv::Mat_<float>& dst, int direct /*= 1*/, bool fillBlock)
        {
            if (src.channels() != 1 )
            {
                return;
            }

            if (direct == 0)
            {
                dst.create(1, src.cols);
            }
            else if (direct == 1)
            {
                dst.create(src.rows, 1);
            }

            if (src.depth() == CV_8U)
            {
                _reduceAvg<uchar>(src, dst, direct, fillBlock);
            }
            else if (src.depth() == CV_32F)
            {
                _reduceAvg<float>(src, dst, direct, fillBlock);
            }
        }

        bool intersect(const K::Line_2& line1, const K::Line_2& line2, K::Point_2& point)
        {
            auto result = CGAL::intersection(line1, line2);

            return CGAL::assign(point, result);
        }

        bool intersect(const K::Line_2& line1, const K::Segment_2& segt2, K::Point_2& point)
        {

            auto result = CGAL::intersection(line1, segt2);

            return CGAL::assign(point, result);
        }

        bool intersect(QLineF& line, QList<QPointF>& segments, QPointF& point)
        {
            K::Line_2 l(as(line.p1()), as(line.p2()));

            for(auto iter = segments.begin() + 1; iter != segments.end(); iter++)
            {
                auto pt1 = as(*(iter - 1));
                auto pt2 = as(*iter);
                
                K::Point_2 pt3;

                if(intersect(l, K::Segment_2(pt1, pt2), pt3))
                {
                    point = to(pt3);
                    return true;
                }
            }

            return false;
        }

        bool intersect(QLineF& line1, QLineF& line2, QPointF& point)
        {
            K::Line_2 l1(as(line1.p1()), as(line1.p2()));
            K::Line_2 l2(as(line2.p1()), as(line2.p2()));

            K::Point_2 pt3;

            if(intersect(l1,l2, pt3))
            {
                point = to(pt3);
                return true;
            }

            return false;
        }

        double  y_at_profile(const QList<QPointF> profile, const double x)
        {
            int closeIndex = 0;

            if (profile.size() == 0)
            {
                return -1;
            }

            while (profile[closeIndex].x() > x)
            {
                closeIndex--;

                if (closeIndex < 0)
                {
                    return -1;
                }
            }

            while (profile[closeIndex + 1].x() < x)
            {
                closeIndex++;

                if (closeIndex >= profile.size() - 1)
                {
                    return -1;
                }
            }

            Geom::Line2d line(Geom::Point2d(profile[closeIndex].x(), profile[closeIndex].y()), Geom::Point2d(profile[closeIndex + 1].x(), profile[closeIndex + 1].y()));

            double y = line.y_at_x(x);
            if (!isnan(y))
            {
                return y;
            }
            else
            {
                return -1;
            }
        }

        float  medianX(const QList<QPointF>& points)
        {
            if (points.empty())
            {
                return 0;
            }

            QList<QPointF> sample;

            sample = points;

            std::sort(std::begin(sample), std::end(sample), [](const QPointF& p1, const QPointF& p2) { return p1.x() < p2.x(); });

            return sample[sample.size() / 2].x();
        }

        float  medianY(const QList<QPointF>& points)
        {
            if (points.empty())
            {
                return 0;
            }

            QList<QPointF> sample;

            sample = points;

            std::sort(std::begin(sample), std::end(sample), [](const QPointF& p1, const QPointF& p2) { return p1.y() < p2.y(); });

            return sample[sample.size() / 2].y();
        }

        float  medianV(const QList<double>& values)
        {
            if (values.empty())
            {
                return 0;
            }

            QList<double> sample;

            sample = values;

            std::sort(std::begin(sample), std::end(sample), [](const double& p1, const double& p2) { return p1 < p2; });

            return sample[sample.size() / 2];
        }

        void calCPCPK(QList<double> values, double& cp, double& cpk, const double tolMax, const double tolMin)
        {
            cp = 0.0;
            cpk = 0.0;
            if (values.size() < 2)
            {
                return;
            }

            double avg = std::accumulate(std::begin(values), std::end(values), 0.0) / values.size();
            double sigma = sqrt(std::accumulate(std::begin(values), std::end(values), 0.0, [=](double v1, double v2) { return v1 + pow(v2 - avg, 2); }) / values.size());

            if (sigma > 0)
            {
                cp = (tolMax - tolMin) / (6 * sigma);
                cpk = qMin<double>((tolMax - avg), (avg - tolMin)) / (3 * sigma);
            }
        }

        void calFPY(QList<double> values, double& okRatio, const double tolMax, const double tolMin)
        {
            okRatio = 0.0;

            if (values.size() > 0)
            {
                okRatio = 1.0 * std::accumulate(std::begin(values), std::end(values), 0, [=](int v1, double v2) { if (v2 >= tolMin && v2 <= tolMax) { return v1 + 1; } else { return v1; } }) / values.size();
            }
        }

        void filter(QList<QPointF>& profile, bool sideCorrect, double thresh,  bool lpf, int filterLen)
        {
           
            for (auto iter = profile.begin(); iter != profile.end(); )
            {
                if (iter->y() > 0.1 && iter->y() < 18.0)
                {
                    break;
                }
                else
                {
                    iter = profile.erase(iter);
                } 
            }

            for (auto iter = profile.rbegin(); iter != profile.rend(); )
            {
                if (iter->y()  > 0.1 && iter->y() < 18)
                {
                    break;
                }
                else
                {
                    iter = QList<QPointF>::reverse_iterator(profile.erase((++iter).base()));
                }
            }

            if (sideCorrect && profile.size() > 20)
            {
                for (auto iter = profile.begin(); iter != profile.end(); )
                {
                    if (iter->y() <= 0.1 || iter->y() > 18.0)
                    {
                        iter = profile.erase(iter);
                        continue;
                    }
                  
                    Base::List<Geom::Point2d> begPoints;
                    Geom::Line2d begSideLine(CGAL::ORIGIN, CGAL::ORIGIN);

                    for (auto iter2 = iter; iter2 != profile.end(); iter2++)
                    {
                        if (begPoints.size() >= 10)
                        {
                            break;
                        }

                        begPoints.push_back(Geom::Point2d(iter2->x(), iter2->y()));
                    }

                    Proc::FitLineD::execute(begPoints, begSideLine);

                    if (!begSideLine.is_degenerate())
                    {
                        double error = 0.0;
                        for (int i =0; i< begPoints.size(); i++)
                        {
                            error += sqrt(CGAL::squared_distance(begSideLine, begPoints[i]));
                        }

                        if (error > thresh)
                        {
                            iter = profile.erase(iter);
                        }
                        else
                        {
                            break;
                        }
                    }
                    else
                    {
                        iter = profile.erase(iter);
                    }

                }

                for (auto iter = profile.rbegin(); iter != profile.rend(); )
                {
                    if (iter->y() <= 0.1)
                    {
                        iter = QList<QPointF>::reverse_iterator(profile.erase((++iter).base()));
                        continue;
                    }
                    Base::List<Geom::Point2d> endPoints;
                    Geom::Line2d endSideLine(CGAL::ORIGIN, CGAL::ORIGIN);

                    for (auto iter2 = iter; iter2 != profile.rend(); iter2++)
                    {
                        if (endPoints.size() >= 10)
                        {
                            break;
                        }

                        endPoints.push_back(Geom::Point2d(iter2->x(), iter2->y()));
                    }

                    Proc::FitLineD::execute(endPoints, endSideLine);

                    if (!endSideLine.is_degenerate())
                    {
                        double error = 0.0;
                        for (int i = 0; i < endPoints.size(); i++)
                        {
                            error += sqrt(CGAL::squared_distance(endSideLine, endPoints[i]));
                        }

                        if (error > thresh)
                        {
                            iter = QList<QPointF>::reverse_iterator(profile.erase((++iter).base()));
                        }
                        else
                        {
                            break;
                        }
                    }
                    else
                    {
                        iter = QList<QPointF>::reverse_iterator(profile.erase((++iter).base()));
                    }
               }
           }
        
            if (lpf)
            {
                int ps = profile.size();

                if (ps > filterLen * 2)
                {
                    for (int i = 0; i < ps; i++)
                    {
                        QList<double> vlist;
                        for (int j = qMax(i - filterLen, 0); j <= qMin(i + filterLen, ps - 1); j++)
                        {
                            vlist.push_back(profile[j].y());
                        }

//                        qSort(vlist);
                        std::sort(vlist.begin(), vlist.end());
                        profile[i].setY(vlist[vlist.size() / 2]);
                    }
                }
                
                /*
               
                QVector<float> pys;

                for (auto& p : profile)
                {
                    pys.append(p.y());
                }

                cv::Mat_<float> line = cv::Mat_<float>(1, pys.size(), pys.data());
                cv::Mat_<float> amp, phase;
               // cv::blur(mat, mat, cv::Size(11, 11));

                QList<int> cComositeIndex;
                if (true)
                {
                    for (int i = 1; i <= qMin(filterLen, line.cols); i++)
                    {
                        cComositeIndex << i;
                    }
                }

                if (DFT_1D(line, amp, phase))
                {
                    float* lineD = (float*)line.data;
                    float* ampD = (float*)amp.data;
                    float* phaseD = (float*)phase.data;

                    double maxValue = -100.0;
                    double maxIndex = 0.0;
                    double min = 0.0;
                    double max = 0.0;

                    for (int i = 0; i < amp.cols; i++)
                    {
                        float sum = 0.0;
                        for (int j : cComositeIndex)
                        {
                            if (ampD[j] > 0.001)
                            {
                                sum += ampD[j] * cos(2 * M_PI*j*i / amp.cols + phaseD[j]);
                            }
                        }

                        lineD[i] = sum + ampD[0] / 2 * cos(phaseD[0]);
                    }
                }
                for (int i = 0; i < pys.size(); i++)
                {
                    auto v = pys[i];
                    profile[i].setY(v);
                }
                */
            }
        }

        ///////////////////////////////
        template <typename T>
        void reduce(Geom::Image& src, Geom::Image& dst)
        {
            dst = Geom::Image(1, src.cols, src.depth());
            T* srcData = (T*)src.data;
            T* dstData = (T*)dst.data;
            T* data;
            int colCount = src.cols;
            int rowCount = src.rows;
            T sum = 0;
            int count = 0;
            for (int c = 0; c < colCount; c++)
            {
                dstData[c] = -1;
                data = srcData + c;
                sum = 0;
                count = 0;
                for (int r = 0; r < rowCount; r++)
                {
                    if (*data >= -1000)
                    {
                        sum += *data;
                        count++;
                    }
                    data += colCount;
                }
                if (count > 0)
                {
                    dstData[c] = sum / count;
                }
            }

            for (int i = 0; i < colCount; i++)
            {
                if (dstData[i] >= 0)
                {
                    for (int j = 0; j < i; j++)
                    {
                        dstData[j] = dstData[i];
                    }
                    break;
                }
            }

            for (int i = colCount - 1; i >= 0; i--)
            {
                if (dstData[i] >= 0)
                {
                    for (int j = colCount - 1; j > i; j--)
                    {
                        dstData[j] = dstData[i];
                    }
                    break;
                }
            }
        }

        cv::Mat getFilterKernel(int extend)
        {
            cv::Mat kernel = cv::Mat::zeros(1, extend * 2 + 1, CV_8S);
            char* data = (char*)kernel.data;
            for (int i = 0; i < extend; i++)
            {
                *(data++) = -1;
            }
            data++;
            for (int i = 0; i < extend; i++)
            {
                *(data++) = +1;
            }
            return kernel;
        }

        int executeU08(const Geom::ImageT<uchar>& source, const Geom::Segment2f& profile, int sampleWidth, double sigma, int minAmp, int maxAmp, int extend, int method, int interp, int transition, FindPointsF::ResultType& points)
        {
            auto p1 = profile.point(0);
            auto p2 = profile.point(1);

            // calculate the length of the line
            auto l = cvRound(sqrt(CGAL::squared_distance(p1, p2)));

            // calculate the angle of the line
            float a = atan2(p2.y() - p1.y(), p2.x() - p1.x());

            int s = -1;
            int e = -1;

            float cosa = cos(a);
            float sina = sin(a);

            float x = cvRound(p1.x());
            float y = cvRound(p1.y());

            // sample pixels along the line perpendicular to the profile line
            Geom::Image samples(2 * sampleWidth - 1, l, CV_16S);

            for (int j = 0; j < l; j++)
            {
                for (int i = -sampleWidth + 1; i < sampleWidth; i++)
                {
                    float x2 = x + cosa*j - sina*i;
                    float y2 = y + sina*j + cosa*i;
                    int c = cvFloor(x2);
                    int r = cvFloor(y2);
                    if (c < 0 || c >= source.cols || r < 0 || r >= source.rows)
                    {
                        samples.at<short>(i + sampleWidth - 1, j) = -1;
                        continue;
                    }
                    switch (interp)
                    {
                        // nearest neighbor pixel interpolation
                    case 0:
                        samples.at<short>(i + sampleWidth - 1, j) = source.at<uchar>(r, c);
                        break;
                        // bilinear interpolation
                    case 1:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<short>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<short>(i + sampleWidth - 1, j) = Algo::bilinear<uchar, uchar>(source, x2, y2);
                        break;
                        // bicubic interpolation, NOT IMPLEMENT YET
                    case 2:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<short>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<short>(i + sampleWidth - 1, j) = Algo::cubic<uchar, uchar>(source, x2, y2);
                        break;
                    }
                    if (s == -1)
                    {
                        s = j;
                    }
                    else
                    {
                        e = j;
                    }
                }
            }

            if (s == -1 || e == -1)
            {
                return -1;
            }

            // calculate the mean of sample pixels as pixel on the profile line
            Geom::Image AV;
            // cv::reduce(samples, AV, 0, CV_REDUCE_AVG, CV_16S);
            reduce<short>(samples, AV);

            if (sigma != 1.0)
            {
                cv::GaussianBlur(AV, AV, cv::Size(), sigma); // apply gaussian smoothing
            }

            if (extend > 1)
            {
                cv::filter2D(AV, AV, AV.depth(), getFilterKernel(extend));
            }

            // calculate 1st order gradient and 2nd order gradient along the profile line

            Geom::Image G1;
            cv::Sobel(AV, G1, CV_16S, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(AV, G1, CV_16S, 1, 0, 1, 1);
            //cv::Sobel(AV, G1, CV_16S, 1, 0, 1, 0.5);
            //cv::Scharr(AV, G1, CV_16S, 1, 0, 0.5);

            if (G1.empty())
            {
                return -1;
            }

            Geom::Image G2;
            cv::Sobel(G1, G2, CV_16S, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(G1, G2, CV_16S, 1, 0, 1, 1);
            //cv::Sobel(G1, G2, CV_16S, 1, 0, 1, 0.5);
            //cv::Scharr(G1, G2, CV_16S, 1, 0, 0.5);

            typedef std::tuple<float, int, int> item;

            // locate the candidate measuring positions by finding peaks in 2nd order gradient
            Base::List<item> cols;

            for (int i = 1; i < G2.cols - 1; ++i)
            {
                if (transition == 0 || transition == 1)
                {
                    int k = i;
                    // test if the pixel is transiting from light to dark
                    if (G2.at<short>(0, i) < 0)
                    {
                        int j = i + 1;
                        // skip the flat area
                        while (i < (G2.cols - 1) && G2.at<short>(0, ++i) == 0)
                        {
                        }
                        // test if the pixel is transiting from dark to light
                        // this is a light peak that possibly a mesuring position if true
                        if (G2.at<short>(0, i) > 0 || i >= (G2.cols - 1))
                        {
                            int amp = std::min(abs(G1.at<short>(0, j)), 255);
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                short k1v = abs(G1.at<short>(0, k - 1));
                                short k2v = abs(G1.at<short>(0, k));
                                short k3v = (k + 1 < G1.cols) ? abs(G1.at<short>(0, k + 1)) : 0;
                                short k4v = (k + 2 < G1.cols) ? abs(G1.at<short>(0, k + 2)) : 0;
                                float p;

                                if (k2v > k3v)
                                {
                                    if (2 * k2v - k1v - k3v == 0)
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (double)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (k2v < k3v)
                                {
                                    if (2 * k3v - k2v - k4v == 0)
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (double)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5 + k;
                                }

                                cols.push_back(std::make_tuple<float, int, int>(std::move(p), 0, std::move(amp)));
                            }
                        }
                        i--;
                    }
                }
                if (transition == 0 || transition == 2)
                {
                    int k = i;
                    // test if the pixel is transiting from dark to light
                    if (G2.at<short>(0, i) > 0)
                    {
                        int j = i;
                        // skip the flat area
                        while (i < (G2.cols - 1) && G2.at<short>(0, ++i) == 0)
                        {
                        }
                        // test if the pixel is transiting from light to dark
                        // this is a dark peak that possibly a mesuring position if true
                        if (G2.at<short>(0, i) < 0 || i >= (G2.cols - 1))
                        {
                            int amp = std::min(abs(G1.at<short>(0, j)), 255);
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                short k1v = abs(G1.at<short>(0, k - 1));
                                short k2v = abs(G1.at<short>(0, k));
                                short k3v = (k + 1 < G1.cols) ? abs(G1.at<short>(0, k + 1)) : 0;
                                short k4v = (k + 2 < G1.cols) ? abs(G1.at<short>(0, k + 2)) : 0;
                                float p;

                                if (k2v > k3v)
                                {
                                    if (2 * k2v - k1v - k3v == 0)
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (double)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (k2v < k3v)
                                {
                                    if (2 * k3v - k2v - k4v == 0)
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (double)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5 + k;
                                }

                                cols.push_back(std::make_tuple<float, int, int>(std::move(p), 1, std::move(amp)));
                            }
                        }
                        i--;
                    }
                }
            }

            if (cols.empty())
            {
                return -1;
            }

            item col;

            switch (method)
            {
            case 0:
                col = cols.front();
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 1:
                col = cols.back();
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 2:
                col = *std::max_element(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 3:
                for (int i = 0; i < cols.size(); i++)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 4:
                for (int i = cols.size() - 1; i >= 0; i--)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 5:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) > std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            case 6:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa * std::get<0>(col), y + sina * std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            }

            return 0;
        }

        int executeU08(const Geom::ImageT<uchar>& source, const Geom::Segment2d& profile, int sampleWidth, double sigma, int minAmp, int maxAmp, int extend, int method, int interp, int transition, FindPointsD::ResultType& points)
        {
            auto p1 = profile.point(0);
            auto p2 = profile.point(1);

            // calculate the length of the line
            auto l = cvRound(sqrt(CGAL::squared_distance(p1, p2)));

            // calculate the angle of the line
            double a = atan2(p2.y() - p1.y(), p2.x() - p1.x());

            double x = cvRound(p1.x());
            double y = cvRound(p1.y());

            int s = -1;
            int e = -1;

            double cosa = cos(a);
            double sina = sin(a);

            // sample pixels along the line perpendicular to the profile line
            Geom::Image samples(2 * sampleWidth - 1, l, CV_16S);
            for (int j = 0; j < l; j++)
            {
                for (int i = -sampleWidth + 1; i < sampleWidth; i++)
                {
                    double x2 = x + cosa*j - sina*i;
                    double y2 = y + sina*j + cosa*i;
                    int c = cvFloor(x2);
                    int r = cvFloor(y2);
                    if (c < 0 || c >= source.cols || r < 0 || r >= source.rows)
                    {
                        samples.at<short>(i + sampleWidth - 1, j) = -1;
                        continue;
                    }
                    switch (interp)
                    {
                        // nearest neighbor pixel interpolation
                    case 0:
                        samples.at<short>(i + sampleWidth - 1, j) = source.at<uchar>(r, c);
                        break;
                        // bilinear interpolation
                    case 1:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<short>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<short>(i + sampleWidth - 1, j) = Algo::bilinear<uchar, uchar>(source, x2, y2);
                        break;
                        // bicubic interpolation, NOT IMPLEMENT YET
                    case 2:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<short>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<short>(i + sampleWidth - 1, j) = Algo::cubic<uchar, uchar>(source, x2, y2);
                        break;
                    }
                    if (s == -1)
                    {
                        s = j;
                    }
                    else
                    {
                        e = j;
                    }
                }
            }

            if (s == -1 || e == -1)
            {
                return -1;
            }

            // calculate the mean of sample pixels as pixel on the profile line
            Geom::Image AV;
            // cv::reduce(samples, AV, 0, CV_REDUCE_AVG, CV_16S);
            reduce<short>(samples, AV);

            if (sigma != 1.0)
            {
                cv::GaussianBlur(AV, AV, cv::Size(), sigma); // apply gaussian smoothing
            }

            if (extend > 1)
            {
                cv::filter2D(AV, AV, AV.depth(), getFilterKernel(extend));
            }

            // calculate 1st order gradient and 2nd order gradient along the profile line

            Geom::Image G1;
            cv::Sobel(AV, G1, CV_16S, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(AV, G1, CV_16S, 1, 0, 1, 1);
            //cv::Sobel(AV, G1, CV_16S, 1, 0, 1, 0.5);
            //cv::Scharr(AV, G1, CV_16S, 1, 0, 0.5);

            if (G1.empty())
            {
                return -1;
            }

            Geom::Image G2;
            cv::Sobel(G1, G2, CV_16S, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(G1, G2, CV_16S, 1, 0, 1, 1);
            //cv::Sobel(G1, G2, CV_16S, 1, 0, 1, 0.5);
            //cv::Scharr(G1, G2, CV_16S, 1, 0, 0.5);

            typedef std::tuple<double, int, int> item;

            // locate the candidate measuring positions by finding peaks in 2nd order gradient
            Base::List<item> cols;

            for (int i = 1; i < G2.cols - 1; ++i)
            {
                if (transition == 0 || transition == 1)
                {
                    int k = i;
                    // test if the pixel is transiting from light to dark
                    if (G2.at<short>(0, i) < 0)
                    {
                        int j = i + 1;
                        // skip the flat area
                        while (i < (G2.cols - 1) && G2.at<short>(0, ++i) == 0)
                        {
                        }
                        // test if the pixel is transiting from dark to light
                        // this is a light peak that possibly a mesuring position if true
                        if (G2.at<short>(0, i) > 0 || i >= (G2.cols - 1))
                        {
                            int amp = std::min(abs(G1.at<short>(0, j)), 255);
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                short k1v = abs(G1.at<short>(0, k - 1));
                                short k2v = abs(G1.at<short>(0, k));
                                short k3v = (k + 1 < G1.cols) ? abs(G1.at<short>(0, k + 1)) : 0;
                                short k4v = (k + 2 < G1.cols) ? abs(G1.at<short>(0, k + 2)) : 0;
                                double p;

                                if (k2v > k3v)
                                {
                                    if (2 * k2v - k1v - k3v == 0)
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (double)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (k2v < k3v)
                                {
                                    if (2 * k3v - k2v - k4v == 0)
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (double)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5 + k;
                                }

                                cols.push_back(std::make_tuple<double, int, int>(std::move(p), 0, std::move(amp)));
                            }
                        }
                        i--;
                    }
                }
                if (transition == 0 || transition == 2)
                {
                    int k = i;
                    // test if the pixel is transiting from dark to light
                    if (G2.at<short>(0, i) > 0)
                    {
                        int j = i;
                        // skip the flat area
                        while (i < (G2.cols - 1) && G2.at<short>(0, ++i) == 0)
                        {
                        }
                        // test if the pixel is transiting from light to dark
                        // this is a dark peak that possibly a mesuring position if true
                        if (G2.at<short>(0, i) < 0 || i >= (G2.cols - 1))
                        {
                            int amp = std::min(abs(G1.at<short>(0, j)), 255);
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                short k1v = abs(G1.at<short>(0, k - 1));
                                short k2v = abs(G1.at<short>(0, k));
                                short k3v = (k + 1 < G1.cols) ? abs(G1.at<short>(0, k + 1)) : 0;
                                short k4v = (k + 2 < G1.cols) ? abs(G1.at<short>(0, k + 2)) : 0;
                                double p;

                                if (k2v > k3v)
                                {
                                    if (2 * k2v - k1v - k3v == 0)
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (double)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (k2v < k3v)
                                {
                                    if (2 * k3v - k2v - k4v == 0)
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (double)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5 + k;
                                }

                                cols.push_back(std::make_tuple<double, int, int>(std::move(p), 1, std::move(amp)));
                            }
                        }
                        i--;
                    }
                }
            }

            if (cols.empty())
            {
                return -1;
            }

            item col;

            switch (method)
            {
            case 0:
                col = cols.front();
                points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 1:
                col = cols.back();
                points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 2:
                col = *std::max_element(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 3:
                for (int i = 0; i < cols.size(); i++)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 4:
                for (int i = cols.size() - 1; i >= 0; i--)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 5:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) > std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            case 6:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa * std::get<0>(col), y + sina * std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            }

            return 0;
        }
        /////////////////////
        int FindLow32(const Geom::ImageT<float>& source, const Geom::Segment2f& profile, int sampleWidth, double sigma, int minAmp, int maxAmp, int extend, int method, int interp, FindPointsF::ResultType& points)
        {
            auto p1 = profile.point(0);
            auto p2 = profile.point(1);

            // calculate the length of the line
            auto l = cvRound(sqrt(CGAL::squared_distance(p1, p2)));

            // calculate the angle of the line
            float a = atan2(p2.y() - p1.y(), p2.x() - p1.x());

            int s = -1;
            int e = -1;

            float cosa = cos(a);
            float sina = sin(a);

            float x = cvRound(p1.x());
            float y = cvRound(p1.y());


            // sample pixels along the line perpendicular to the profile line
            Geom::ImageT<float> samples(2 * sampleWidth - 1, l);
            for (int j = 0; j < l; j++)
            {
                for (int i = -sampleWidth + 1; i < sampleWidth; i++)
                {
                    float x2 = x + cosa*j - sina*i;
                    float y2 = y + sina*j + cosa*i;
                    int c = cvFloor(x2);
                    int r = cvFloor(y2);
                    if (c < 0 || c >= source.cols || r < 0 || r >= source.rows)
                    {
                        samples.at<float>(i + sampleWidth - 1, j) = -1;
                        continue;
                    }
                    float a;
                    switch (interp)
                    {
                        // nearest neighbor pixel interpolation
                    case 0:
                        samples.at<float>(i + sampleWidth - 1, j) = source.at<float>(r, c);
                        break;
                        // bilinear interpolation
                    case 1:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<float>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        a = Algo::bilinear<float, float>(source, x2, y2);
                        samples.at<float>(i + sampleWidth - 1, j) = Algo::bilinear<float, float>(source, x2, y2);
                        break;
                        // bicubic interpolation, NOT IMPLEMENT YET
                    case 2:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<float>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<float>(i + sampleWidth - 1, j) = Algo::cubic<float, float>(source, x2, y2);
                        break;
                    }
                    if (s == -1)
                    {
                        s = j;
                    }
                    else
                    {
                        e = j;
                    }
                }
            }

            if (s == -1 || e == -1)
            {
                return -1;
            }

            // calculate the mean of sample pixels as pixel on the profile line
            Geom::Image AV;

            reduce<float>(samples, AV);

            if (sigma != 1.0)
            {
                cv::GaussianBlur(AV, AV, cv::Size(), 2.0); // apply gaussian smoothing
            }

            if (extend > 1)
            {
                cv::filter2D(AV, AV, AV.depth(), getFilterKernel(extend));
            }

            // calculate 1st order gradient and 2nd order gradient along the profile line

            Geom::Image G1;
            cv::Sobel(AV, G1, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);


            if (G1.empty())
            {
                return -1;
            }

            Geom::Image G2;
            cv::Sobel(G1, G2, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);

            typedef std::tuple<int, int, int> item;

            // locate the candidate measuring positions by finding peaks in 2nd order gradient
            Base::List<item> cols;
            for (int i = 0; i < G1.cols ; ++i)
            {
                if (G1.at<float>(0, i) >= 0)
                {
                    continue;
                }
                else {
                    float aut = fabs(G1.at<float>(0, i));
                    if (aut > minAmp && aut < maxAmp)
                    {
                        int k = i;
                        while (G1.at<float>(0, i) < 0)
                        {
                            if (i < G1.cols )
                            {
                                if (i > 0)
                                {
                                    if (G2.at<float>(0, i - 1) * G2.at<float>(0, i ) <= 0)
                                    {
                                        k = i;
                                    }
                                }
                                else {
                                    if (G2.at<float>(0, i ) * G2.at<float>(0, i + 1) <= 0)
                                    {
                                        k = i;
                                    }
                                }
                                if (i < G1.cols - 1)
                                {
                                    i++;
                                }
                                else {
                                    return 0;
                                }
                            }
                            else {
                                return 0;
                            }

                        }

                        int k1 = i ;
                        if (i < G1.cols - 1)
                        {
                            i++;
                        }
                        if (k1 < k)
                        {
                            continue;
                        }
                        int k2 = i;
                        while (G1.at<float>(0, i) >= 0)
                        {
                            if (i < G1.cols)
                            {
                                if (i > 0)
                                {
                                    if (G2.at<float>(0, i - 1) * G2.at<float>(0, i ) <= 0)
                                    {
                                        k2 = i-1 ;
                                        break;
                                    }
                                }
                                else {
                                    if (G2.at<float>(0, i ) * G2.at<float>(0, i + 1) <= 0)
                                    {
                                        k2 = i;
                                    }
                                }
                                if (i < G1.cols - 1)
                                {
                                    i++;
                                }
                                else {
                                    break;
                                }
                            }
                            else {
                                return 0;
                            }
                        }
                        if (k1 > k2)
                        {
                            continue;
                        }
                        if (k1 == k2 && k1 == G1.cols - 1)
                        {
                            continue;
                        }
                        cols.push_back(std::make_tuple(k, k1, k2));
                    }
                }
            }

            if (cols.empty())
            {
                return -1;
            }

            item col;

            switch (method)
            {
            case 0:
                col = cols.front();
                for (auto iter = cols.begin(); iter != cols.end(); ++iter)
                {
                    if (AV.at<float>(0, abs(std::get<1>(*iter))) < AV.at<float>(0, abs(std::get<1>(col))))
                    {
                        col = *iter;
                    }
                }
                if (abs(std::get<2>(col) - std::get<0>(col)) < G1.cols/2 && abs(std::get<2>(col) - std::get<0>(col)) > G1.cols / 20)
                {
                    if (abs(std::get<2>(col) - std::get<0>(col)) >= 2 && abs(AV.at<float>(0, abs(std::get<2>(col))) - AV.at<float>(0, abs(std::get<1>(col)))) > abs(AV.at<float>(0, abs(std::get<0>(col))) - AV.at<float>(0, abs(std::get<1>(col)))) ? abs(AV.at<float>(0, abs(std::get<2>(col))) - AV.at<float>(0, abs(std::get<1>(col)))) : abs(AV.at<float>(0, abs(std::get<0>(col))) - AV.at<float>(0, abs(std::get<1>(col)))) > 0.09)
                    {
                        double p;
                        p = std::get<1>(col) - 1 + (double)1.3 * abs(G1.at<float>(0, std::get<1>(col) - 1) ) / (abs(G1.at<float>(0, std::get<1>(col) - 1)) + abs(G1.at<float>(0, std::get<1>(col)) ));
                        points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa * p, y + sina * p), 0, 0));
                    }
                    else {
                        points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa * std::get<1>(col), y + sina * std::get<1>(col)), 0, 0));
                    }
                }
                break;
            case 1:
                col = cols.back();
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<1>(col), y + sina*std::get<1>(col)), 0, 0));
                break;
            }

            return 0;
        }
        ///////////////////////////

        int executeF32(const Geom::ImageT<float>& source, const Geom::Segment2f& profile, int sampleWidth, double sigma, int minAmp, int maxAmp, int extend, int method, int interp, int transition, FindPointsF::ResultType& points)
        {
            auto p1 = profile.point(0);
            auto p2 = profile.point(1);

            // calculate the length of the line
            auto l = cvRound(sqrt(CGAL::squared_distance(p1, p2)));

            // calculate the angle of the line
            float a = atan2(p2.y() - p1.y(), p2.x() - p1.x());

            int s = -1;
            int e = -1;

            float cosa = cos(a);
            float sina = sin(a);

            float x = cvRound(p1.x());
            float y = cvRound(p1.y());

            
            // sample pixels along the line perpendicular to the profile line
            Geom::ImageT<float> samples(2 * sampleWidth - 1, l);
            for (int j = 0; j < l; j++)
            {
                for (int i = -sampleWidth + 1; i < sampleWidth; i++)
                {
                    float x2 = x + cosa*j - sina*i;
                    float y2 = y + sina*j + cosa*i;
                    int c = cvFloor(x2);
                    int r = cvFloor(y2);
                    if (c < 0 || c >= source.cols || r < 0 || r >= source.rows)
                    {
                        samples.at<float>(i + sampleWidth - 1, j) = -1;
                        continue;
                    }
                    float a;
                    switch (interp)
                    {
                        // nearest neighbor pixel interpolation
                    case 0:
                        samples.at<float>(i + sampleWidth - 1, j) = source.at<float>(r, c);
                        break;
                        // bilinear interpolation
                    case 1:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<float>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        a= Algo::bilinear<float, float>(source, x2, y2);
                        samples.at<float>(i + sampleWidth - 1, j) = Algo::bilinear<float, float>(source, x2, y2);
                        break;
                        // bicubic interpolation, NOT IMPLEMENT YET
                    case 2:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<float>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<float>(i + sampleWidth - 1, j) = Algo::cubic<float, float>(source, x2, y2);
                        break;
                    }
                    if (s == -1)
                    {
                        s = j;
                    }
                    else
                    {
                        e = j;
                    }
                }
            }

            if (s == -1 || e == -1)
            {
                return -1;
            }

            // calculate the mean of sample pixels as pixel on the profile line
            Geom::Image AV;
            // cv::reduce(samples, AV, 0, CV_REDUCE_AVG, CV_32F);
            reduce<float>(samples, AV);

            if (sigma != 1.0)
            {
                cv::GaussianBlur(AV, AV, cv::Size(), sigma); // apply gaussian smoothing
            }

            if (extend > 1)
            {
                cv::filter2D(AV, AV, AV.depth(), getFilterKernel(extend));
            }

            // calculate 1st order gradient and 2nd order gradient along the profile line

            Geom::Image G1;
            cv::Sobel(AV, G1, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(AV, G1, CV_32F, 1, 0, 1, 1);
            //cv::Sobel(AV, G1, CV_32F, 1, 0, 1, 0.5);
            //cv::Scharr(AV, G1, CV_32F, 1, 0, 0.5);

            if (G1.empty())
            {
                return -1;
            }

            Geom::Image G2;
            cv::Sobel(G1, G2, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(G1, G2, CV_32F, 1, 0, 1, 1);
            //cv::Sobel(G1, G2, CV_32F, 1, 0, 1, 0.5);
            //cv::Scharr(G1, G2, CV_32F, 1, 0, 0.5);

            typedef std::tuple<float, int, float> item;

            /*std::ofstream out2("C:/2.txt",std::ios::app);
            if (!out2.is_open())
            {
                return 1;
            }

            for (int i = 0; i < G2.cols; ++i)
            {
                out2 << AV.at<float>(0, i)<<"      "<< abs(G1.at<float>(0, i)) << std::endl;
            }*/
            // locate the candidate measuring positions by finding peaks in 2nd order gradient
            Base::List<item> cols;

            for (int i = 1; i < G2.cols - 1; ++i)
            {
                if (transition == 0 || transition == 1)
                {
                    int k = i;
                    // test if the pixel is transiting from light to dark
                    if (G2.at<float>(0, i) < 0)
                    {
                        int j = i + 1;
                        // skip the flat area
                        while (i < (G2.cols - 1) && qFuzzyIsNull(G2.at<float>(0, ++i)))
                        {
                        }
                        // test if the pixel is transiting from dark to light
                        // this is a light peak that possibly a mesuring position if true
                        if (G2.at<float>(0, i) > 0)
                        {
                            float amp = fabs(G1.at<float>(0, j));
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                float k1v = (G1.at<float>(0, k - 1));
                                float k2v = (G1.at<float>(0, k));
                                float k3v = (k + 1 < G1.cols) ? (G1.at<float>(0, k + 1)) : 0;
                                float k4v = (k + 2 < G1.cols) ? (G1.at<float>(0, k + 2)) : 0;
                                float p;

                                if (fabs(k2v) > fabs(k3v))
                                {
                                    if (qFuzzyIsNull(2 * k2v - k1v - k3v))
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (float)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (fabs(k2v) < fabs(k3v))
                                {
                                    if (qFuzzyIsNull(2 * k3v - k2v - k4v) || qFuzzyIsNull(k4v))
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (float)((-k2v)*(2 * k - 3) + (2 * k3v)*(2 * k - 1) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5f + k;
                                }
                                if (p < -1 || p>l + 1)
                                {
                                    //out2 << "p:" << k + 0.5f<<std::endl;
                                    //out2 << "k:" << k << "  k1v:" << k1v << "  k2v:" << k2v << "  k3v:" << k3v << "  k4v:" << k4v << "  p:" << p <<"  klv:"<< (2 * k3v - k2v - k4v)<< std::endl;
                                    cols.push_back(std::make_tuple<float, int, float>(std::move(k + 0.5f), 0, std::move(amp)));
                                }
                                else {
                                    //out2 << "p:" << p << std::endl;
                                    cols.push_back(std::make_tuple<float, int, float>(std::move(p), 0, std::move(amp)));
                                }

                            }
                        }
                        else {
                            if (method == 5 || method == 6 || method == 2)
                            {
                                cols.push_back(std::make_tuple<float, int, float>(std::move(i), 0, std::move(fabs(G1.at<float>(0, i)))));
                            }

                        }
                        i--;
                    }
                }
                if (transition == 0 || transition == 2)
                {
                    int k = i;
                    // test if the pixel is transiting from dark to light
                    if (G2.at<float>(0, i) > 0)
                    {
                        int j = i;
                        // skip the flat area
                        while (i < (G2.cols - 1) && qFuzzyIsNull(G2.at<float>(0, ++i)))
                        {
                        }
                        // test if the pixel is transiting from light to dark
                        // this is a dark peak that possibly a mesuring position if true
                        if (G2.at<float>(0, i) < 0)
                        {
                            float amp = fabs(G1.at<float>(0, j));
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                float k1v = G1.at<float>(0, k - 1);
                                float k2v = G1.at<float>(0, k);
                                float k3v = (k + 1 < G1.cols) ? G1.at<float>(0, k + 1) : 0;
                                float k4v = (k + 2 < G1.cols) ? G1.at<float>(0, k + 2) : 0;
                                float p;

                                if (fabs(k2v) > fabs(k3v))
                                {
                                    if (qFuzzyIsNull(2 * k2v - k1v - k3v))
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (float)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (fabs(k2v) < fabs(k3v))
                                {
                                    if (qFuzzyIsNull(2 * k3v - k2v - k4v) || qFuzzyIsNull(k4v))
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (float)((-k2v)*(2 * k - 3) + (2 * k3v)*(2 * k - 1) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v); //(float)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5f + k;
                                }

                                if (p < -1 || p>l + 1)
                                {
                                    //out2 << "p:" << k + 0.5f << std::endl;
                                    //out2 << "k:" << k << "  k1v:" << k1v << "  k2v:" << k2v << "  k3v:" << k3v << "  k4v:" << k4v << "  p:" << p << "  klv:" << (2 * k3v - k2v - k4v) << std::endl;
                                    cols.push_back(std::make_tuple<float, int, float>(std::move(k + 0.5f), 1, std::move(amp)));
                                }
                                else {
                                    //out2 << "p:" << p << std::endl;
                                    cols.push_back(std::make_tuple<float, int, float>(std::move(p), 1, std::move(amp)));
                                }

                            }
                        }
                        else {
                            if (method == 5 || method == 6 || method==2 )
                            {
                                cols.push_back(std::make_tuple<float, int, float>(std::move(i), 1, std::move(fabs(G1.at<float>(0, i)))));
                            }
                        }
                        i--;
                    }
                }
            }

            if (cols.empty())
            {
                return -1;
            }
            Base::List<item>::iterator iter,iter1;
            if (cols.size() >= 2)
            {
                iter = cols.begin(); iter1 = iter + 1;
                while (iter1 != cols.end())
                {
                    if (fabs(std::get<0>(*iter1) - std::get<0>(*iter)) <= 2)
                    {
                        if (fabs(std::get<2>(*iter1) - std::get<2>(*iter)) < 0.015)
                        {
                            iter1 = cols.erase(iter1);
                        }
                        else {
                            iter++;
                            iter1++;
                        }
                    }
                    else {
                        iter++;
                        iter1++;
                    }

                }
            }
            

            item col;

            switch (method)
            {
            case 0:
                col = cols.front();
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 1:
                col = cols.back();
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 2:
                col = *std::max_element(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 3:
                for (int i = 0; i < cols.size(); i++)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 4:
                for (int i = cols.size() - 1; i >= 0; i--)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 5:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) > std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            case 6:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsF::Node(Geom::Point2f(x + cosa * std::get<0>(col), y + sina * std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            }

            return 0;
        }

        int executeF32(const Geom::ImageT<float>& source, const Geom::Segment2d& profile, int sampleWidth, double sigma, int minAmp, int maxAmp, int extend, int method, int interp, int transition, FindPointsD::ResultType& points)
        {
            auto p1 = profile.point(0);
            auto p2 = profile.point(1);

            // calculate the length of the line
            auto l = cvRound(sqrt(CGAL::squared_distance(p1, p2)));

            // calculate the angle of the line
            double a = atan2(p2.y() - p1.y(), p2.x() - p1.x());

            double x = cvRound(p1.x());
            double y = cvRound(p1.y());

            int s = -1;
            int e = -1;

            double cosa = cos(a);
            double sina = sin(a);

            // sample pixels along the line perpendicular to the profile line
            Geom::ImageT<float> samples(2 * sampleWidth - 1, l);
            for (int j = 0; j < l; j++)
            {
                for (int i = -sampleWidth + 1; i < sampleWidth; i++)
                {
                    double x2 = x + cosa*j - sina*i;
                    double y2 = y + sina*j + cosa*i;
                    int c = cvFloor(x2);
                    int r = cvFloor(y2);
                    if (c < 0 || c >= source.cols || r < 0 || r >= source.rows)
                    {
                        samples.at<float>(i + sampleWidth - 1, j) = -1;
                        continue;
                    }
                    switch (interp)
                    {
                        // nearest neighbor pixel interpolation
                    case 0:
                        samples.at<float>(i + sampleWidth - 1, j) = source.at<float>(r, c);
                        break;
                        // bilinear interpolation
                    case 1:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<float>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<float>(i + sampleWidth - 1, j) = Algo::bilinear<float, float>(source, x2, y2);
                        break;
                        // bicubic interpolation, NOT IMPLEMENT YET
                    case 2:
                        if (c + 1 >= source.cols || r + 1 >= source.rows)
                        {
                            samples.at<float>(i + sampleWidth - 1, j) = -1;
                            continue;
                        }
                        samples.at<float>(i + sampleWidth - 1, j) = Algo::cubic<float, float>(source, x2, y2);
                        break;
                    }
                    if (s == -1)
                    {
                        s = j;
                    }
                    else
                    {
                        e = j;
                    }
                }
            }

            if (s == -1 || e == -1)
            {
                return -1;
            }

            // calculate the mean of sample pixels as pixel on the profile line
            Geom::Image AV;
            // cv::reduce(samples, AV, 0, CV_REDUCE_AVG, CV_32F);
            reduce<float>(samples, AV);

            if (sigma != 1.0)
            {
                cv::GaussianBlur(AV, AV, cv::Size(), sigma); // apply gaussian smoothing
            }

            if (extend > 1)
            {
                cv::filter2D(AV, AV, AV.depth(), getFilterKernel(extend));
            }

            // calculate 1st order gradient and 2nd order gradient along the profile line

            Geom::Image G1;
            cv::Sobel(AV, G1, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(AV, G1, CV_32F, 1, 0, 1, 1);
            //cv::Sobel(AV, G1, CV_32F, 1, 0, 1, 0.5);
            //cv::Scharr(AV, G1, CV_32F, 1, 0, 0.5);

            if (G1.empty())
            {
                return -1;
            }

            Geom::Image G2;
            cv::Sobel(G1, G2, CV_32F, 1, 0, 3, 1.0, 0.0, cv::BORDER_REPLICATE);
            //cv::Sobel(G1, G2, CV_32F, 1, 0, 1, 1);
            //cv::Sobel(G1, G2, CV_32F, 1, 0, 1, 0.5);
            //cv::Scharr(G1, G2, CV_32F, 1, 0, 0.5);

            typedef std::tuple<double, int, int> item;

            // locate the candidate measuring positions by finding peaks in 2nd order gradient
            Base::List<item> cols;

            for (int i = 1; i < G2.cols - 1; ++i)
            {
                if (transition == 0 || transition == 1)
                {
                    int k = i;
                    // test if the pixel is transiting from light to dark
                    if (G2.at<float>(0, i) < 0)
                    {
                        int j = i + 1;
                        // skip the flat area
                        while (i < (G2.cols - 1) && qFuzzyIsNull(G2.at<float>(0, ++i)))
                        {
                        }
                        // test if the pixel is transiting from dark to light
                        // this is a light peak that possibly a mesuring position if true
                        if (G2.at<float>(0, i) > 0 || i >= (G2.cols - 1))
                        {
                            float amp = abs(G1.at<float>(0, j));
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                float k1v = (G1.at<float>(0, k - 1));
                                float k2v = (G1.at<float>(0, k));
                                float k3v = (k + 1 < G1.cols) ? (G1.at<float>(0, k + 1)) : 0;
                                float k4v = (k + 2 < G1.cols) ? (G1.at<float>(0, k + 2)) : 0;
                                double p;

                                if (k2v > k3v)
                                {
                                    if (qFuzzyIsNull(2 * k2v - k1v - k3v))
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (double)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (k2v < k3v)
                                {
                                    if (qFuzzyIsNull(2 * k3v - k2v - k4v))
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (double)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5 + k;
                                }

                                cols.push_back(std::make_tuple<double, int, int>(std::move(p), 0, std::move(amp)));
                            }
                        }
                        i--;
                    }
                }
                if (transition == 0 || transition == 2)
                {
                    int k = i;
                    // test if the pixel is transiting from dark to light
                    if (G2.at<float>(0, i) > 0)
                    {
                        int j = i;
                        // skip the flat area
                        while (i < (G2.cols - 1) && qFuzzyIsNull(G2.at<float>(0, ++i)))
                        {
                        }
                        // test if the pixel is transiting from light to dark
                        // this is a dark peak that possibly a mesuring position if true
                        if (G2.at<float>(0, i) < 0 || i >= (G2.cols - 1))
                        {
                            float amp = abs(G1.at<float>(0, j));
                            if (amp >= minAmp && amp <= maxAmp)
                            {
                                float k1v = (G1.at<float>(0, k - 1));
                                float k2v = (G1.at<float>(0, k));
                                float k3v = (k + 1 < G1.cols) ? (G1.at<float>(0, k + 1)) : 0;
                                float k4v = (k + 2 < G1.cols) ? (G1.at<float>(0, k + 2)) : 0;
                                double p;

                                if (k2v > k3v)
                                {
                                    if (qFuzzyIsNull(2 * k2v - k1v - k3v))
                                    {
                                        p = k;
                                    }
                                    else
                                    {
                                        p = (double)((-k1v)*(2 * k + 1) + (4 * k2v*k) + (-k3v)*(2 * k - 1)) / (4 * k2v - 2 * k1v - 2 * k3v);
                                    }
                                }
                                else if (k2v < k3v)
                                {
                                    if (qFuzzyIsNull(2 * k3v - k2v - k4v))
                                    {
                                        p = k + 1;
                                    }
                                    else
                                    {
                                        p = (double)((-k2v)*(2 * k + 3) + (2 * k3v)*(2 * k + 2) + (-k4v)*(2 * k + 1)) / (4 * k3v - 2 * k2v - 2 * k4v);
                                    }
                                }
                                else
                                {
                                    p = 0.5 + k;
                                }

                                cols.push_back(std::make_tuple<double, int, int>(std::move(p), 1, std::move(amp)));
                            }
                        }
                        i--;
                    }
                }
            }

            if (cols.empty())
            {
                return -1;
            }

            item col;

            switch (method)
            {
            case 0:
                col = cols.front();
                points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 1:
                col = cols.back();
                points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 2:
                col = *std::max_element(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                break;
            case 3:
                for (int i = 0; i < cols.size(); i++)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 4:
                for (int i = cols.size() - 1; i >= 0; i--)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(cols[i]), y + sina*std::get<0>(cols[i])), std::get<1>(cols[i]), std::get<2>(cols[i])));
                }
                break;
            case 5:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) > std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa*std::get<0>(col), y + sina*std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            case 6:
                std::sort(cols.begin(), cols.end(), [](const item& item1, const item& item2) { return std::get<2>(item1) < std::get<2>(item2); });
                for (auto& col : cols)
                {
                    points.push_back(FindPointsD::Node(Geom::Point2d(x + cosa * std::get<0>(col), y + sina * std::get<0>(col)), std::get<1>(col), std::get<2>(col)));
                }
                break;
            }

            return 0;
        }
    
    int  FindPointsF::execute(const Geom::Image& source, const Geom::Segment2f& profile, int sampleWidth, double sigma, double minAmp, double maxAmp, int extend, int method, int interp, int transition, ResultType& points)
    {
        if (CV_MAT_DEPTH(source.type()) == CV_8U)
        {
            return executeU08(source, profile, sampleWidth, sigma, minAmp, maxAmp, extend, method, interp, transition, points);
        }
        else
        {
            return executeF32(source, profile, sampleWidth, sigma, minAmp, maxAmp, extend, method, interp, transition, points);
        }
    }

    int FindPointsD::execute(const Geom::Image& source, const Geom::Segment2d& profile, int sampleWidth, double sigma, double minAmp, double maxAmp, int extend, int method, int interp, int transition, ResultType& points)
    {
        if (CV_MAT_DEPTH(source.type()) == CV_8U)
        {
            return executeU08(source, profile, sampleWidth, sigma, minAmp, maxAmp, extend, method, interp, transition, points);
        }
        else
        {
            return executeF32(source, profile, sampleWidth, sigma, minAmp, maxAmp, extend, method, interp, transition, points);
        }
    }


    
} }